### Robustness: Interaction between skin group and racial identification


## Regression with interactions on both race and skin color & triple -------
reg_interactions <- function(dep_var) {
  
  ind_controls <- paste("male", "factor(religion_simple)", 
                        "age_grade_distortion", "score_self_esteem", 
                        "score_parents_support", "score_study", sep = "+")
  
  ses_controls <- paste("scores_poverty", "score_neighborhood_quality", 
                        "factor(occ_father_simple_no_job_no_na)",
                        "factor(occ_mother_simple_no_job_no_na)", sep = "+")
  
  class_fe <- paste("factor(class_id)")
  
  
  basic_reg_sc_nw <- paste0(dep_var, "~", "factor(negro_na)*grades_score + skin_group_2*grades_score")
  basic_reg_triple_nw <- paste0(dep_var, "~", "factor(negro_na)*grades_score*skin_group_2")
  
  list_basic_regs <- c(basic_reg_sc_nw, basic_reg_triple_nw)
  # basic regression
  Form_basic <- lapply(list_basic_regs, formula)
  
  # Individual control
  Form_ind_control <- lapply(list_basic_regs, function(s) formula(paste(s, ind_controls, sep = "+")))
  
  # SES control
  Form_ses_control <- lapply(list_basic_regs, function(s) formula(paste(s, ses_controls, sep = "+")))
  
  # Individual + SES control
  Form_ind_ses_control <- lapply(list_basic_regs, function(s) formula(paste(s, ind_controls, ses_controls, sep = "+")))
  
  # Individual control + class FE
  Form_ind_control_class_fe <- lapply(list_basic_regs, function(s) formula(paste(s, ind_controls, class_fe, sep = "+")))
  
  # SES control + class FE
  Form_ses_control_class_fe <- lapply(list_basic_regs, function(s) formula(paste(s, ses_controls, class_fe, sep = "+")))
  
  # Individual + SES control + class FE
  Form_ind_ses_control_class_fe <- lapply(list_basic_regs, function(s) formula(paste(s, ind_controls, ses_controls, class_fe, sep = "+")))
  
  
  list_specifications <- list(Form_basic, Form_ind_control, Form_ses_control, 
                              Form_ind_ses_control, Form_ind_control_class_fe,
                              Form_ses_control_class_fe, Form_ind_ses_control_class_fe) 
  results <- list()
  results_r <- list()
  results_c <- list()
  
  for(i in 1:length(list_specifications)) {
    results[[i]] <- lapply(list_specifications[[i]], function(x) (lm(x, data = data, subset = (race < 4))))
    results_r[[i]] <- lapply(results[[i]], robust_se, "HC3")
    results_c[[i]] <- lapply(results[[i]], cluster_se, "class_id")
  }
  
  results_final <- list(results, results_r, results_c)  
  return(results_final)
}

results_same <- reg_interactions("ssi_same_negro_na_2_std")
results_all <- reg_interactions("ssi_all_race_2_std")
results_other <- reg_interactions("ssi_other_race_negro_na_2_std")

list_regs <- list(results_same[[1]][[7]][[1]], results_same[[1]][[7]][[2]], 
                  results_other[[1]][[7]][[1]], results_other[[1]][[7]][[2]],
                  results_all[[1]][[7]][[1]], results_all[[1]][[7]][[2]])

list_regs_c <- list(results_same[[3]][[7]][[1]], results_same[[3]][[7]][[2]], 
                    results_other[[3]][[7]][[1]], results_other[[3]][[7]][[2]],
                    results_all[[3]][[7]][[1]], results_all[[3]][[7]][[2]])

n_obs_bi <- sapply(list_regs, get_obs)

r_adj_bi <- sapply(list_regs, get_adj_r)


variables_to_omit_nl <- "(age)|(male)|(religion)|(class_id)|(avg_)|(class)|(school)|(missing)|(EM)|
                      |(sh)|(poverty)|(Constant)|(supply)|
|(score_nei)|(score_sel)|(score_pare)|(score_study)|(occ)"


cov_var_labels_interaction_bi <- c("Nonwhite", "Grades", "Dark Skin", 
                                   "Nonwhite*Grades", "Nonwhite*Dark Skin",
                                   "Dark Skin*Grades",
                                   "Nonwhite*Dark Skin*Grades")

## Table S9
stargazer(list_regs_c,
          column.separate = c(2,2,2),
          covariate.labels = cov_var_labels_interaction_bi,
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("Same race", "Other races", "All races"),
          omit = variables_to_omit_nl, 
          omit.stat = c("f", "rsq","ser"),
          add.lines = list(c("Individual controls", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Classroom FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Observation", n_obs_bi),
                           c("Adjusted R2", r_adj_bi)),
          font.size = "tiny", type = "text")

stargazer(list_regs_c,
          column.separate = c(2,2,2),
          covariate.labels = cov_var_labels_interaction_bi,
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("Same race", "Other races", "All races"),
          omit = variables_to_omit_nl, 
          omit.stat = c("f", "rsq","ser"),
          add.lines = list(c("Individual controls", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Classroom FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Observation", n_obs_bi),
                           c("Adjusted R2", r_adj_bi)),
          font.size = "footnotesize", type = "latex", 
          out = "tables/regs_1_interaction.tex",
          column.sep.width = "-5pt",
          label = "regs_1_interaction",
          title = "Racial identity, skin color, grades, and social status")



